Persistence and stability of generalized ribosome flow models with time-varying transition rates

In this paper some important qualitative dynamical properties of generalized ribosome flow models are studied. Ribosome flow models known from the literature are generalized by allowing an arbitrary directed network structure between compartments, and by assuming general time-varying rate functions corresponding to the transitions. Persistence of the dynamics is shown using the chemical reaction network (CRN) representation of the system where the state variables correspond to ribosome density and the amount of free space in the compartments. The L1 contractivity of solutions is also proved in the case of periodic reaction rates having the same period. Further we prove the stability of different compartmental structures including strongly connected ones with entropy-like logarithmic Lyapunov functions through embedding the model into a weakly reversible CRN with time-varying reaction rates in a reduced state space. Moreover, it is shown that different Lyapunov functions may be assigned to the same model depending on the non-unique factorization of the reaction rates. The results are illustrated through several examples with biological meaning including the classical ribosome flow model on a ring.


Introduction
Compartmental models are used to describe and analyze the transport between different containers, called compartments in various natural and technological systems [1,2]. Compartments can be assigned to tissues or organs in pharmacokinetic models, mass containers in process systems, distinct disease states in epidemiological models, road sections in transportation systems or different habitats in ecological models. The modeled objects (molecules, people, vehicles, etc.) can move between compartments obeying the given constraints such as limits of directions, flow rates, or capacities. A fundamental feature of compartmental models is that each modeled object can be present in exactly one compartment at a given time. Naturally, compartmental models written in the original physical coordinates belong to the class of nonnegative systems for which the nonnegative orthant is invariant with respect to the dynamics [3,4]. This special property supports the dynamical analysis and control design in several ways. The controllability, observability, realizability and identifiability of mainly linear compartmental system are addressed in [5]. An excellent overview of the qualitative dynamical properties of general compartmental systems can be found in [6]. The dynamical modeling of the mRNA translation process has been in the focus of research since the second half of the 20th century (see, e.g. [7][8][9]). The first large scale analysis of gene translation through the so-called ribosome flow model (RFM) was presented in [10], where the applied second order nonnegative and nonlinear model based on the principle of Totally Asymmetric Exclusion [11] was able to capture the most important dynamical features of the translation process. Also in [10], the RFM model was validated through biological data obtained from three different organisms, and it was clearly shown that its predictive power is superior to several other popular techniques. In [12] the RFM was equipped with an appropriate input-output pair, and it was shown that after applying an affine positive output feedback, the system had a unique equilibrium point which is globally stable in the bounded operating domain. A circular RFM structure was analyzed in [13], where the authors proved using the theory of cooperative systems that the system has a continuum of equilibria, but each equilibrium is globally asymptotically stable within the equivalence class of trajectories determined by the initial conditions. The stability of periodic solutions was also shown. In [14] a bounded pool of free ribosomes was added to the RFM generating a competition among the arbitrary number of mRNA molecules for ribosomes. This generates a special network structure for RFM subsystems, for which the uniqueness and stability of equilibria together with the properties of periodic solutions were proved, too. Different compartment sizes of the RFM were assumed in [15], and it was shown that this modification does not change the favorable dynamical properties of the system. In [16], the ribosome flow model with Langmuir kinetics (RFMLK) is introduced, and a network structure is constructed with RFMLK subsystems connected through a pool. Among other results, it is shown that the trajectories of such a network always converge to a unique equilibrium.
Chemical reaction networks (CRNs) also called kinetic systems can be considered as universal descriptors of nonlinear dynamics, especially that of nonnegative systems [17]. Since the 1970's the theory of CRNs has been intensively studied, and there are several fundamental results on the relation between network structure/parametrization and dynamical properties [18]. The stability of mass-action type CRNs is most often analyzed using an entropy-like logarithmic Lyapunov function, originally called a "pseudo-Helmholtz function" in [19]. Probably the most well-known conjecture of chemical reaction network theory is the "Global attractor conjecture" according to which complex balanced kinetic systems are globally stable with respect to the nonnegative orthant with the logarithmic Lyapunov function [20]. This conjecture was proved for complex balanced reaction networks with a reaction graph of one component [21]. One of the most important results from the point of view of this paper is [22] studying zero deficiency networks, where the allowed kinetics is more general than mass action, the rate coefficients can be time-varying, and the logarithmic Lyapunov function is also generalized. The Lyapunov-function-based stability analysis of RFMs is mentioned as an important problem in [23], which will be addressed in this paper using the CRN representation of the system. In [24] a so-called Max-Min type robust Lyapunov function composed of piecewise linear terms was constructed for a tubular RFM with mass action kinetics.
It is interesting to mention that mathematical models which are equivalent to RFMs can also be obtained through a special finite volume spatial discretization of widely used flow models in PDE form [25,26]. These models also have a transparent representation in CRN form supporting further dynamical analysis. An arbitrary directed graph structure of such models with general time-invariant kinetics was considered in [27], where the existence and uniqueness of equlibria, persistence and contractivity (non-expansive property) of the solutions was shown using the theory of Petri nets, compartmental systems, and earlier results on RFMs. The stability of this model class with logarithmic Lyapunov functions was shown in [28], while a port-Hamiltonian description was given in [29].
Based on the above overview, the aim of this paper is to extend the results of [27,28] in the following respects: considering even more general kinetics with explicit time-dependence, the qualitative analysis of periodic solutions, and finally, stability analysis with a family of different logarithmic Lyapunov functions.
The structure of the paper is the following. Section 2 contains the applied mathematical notations for compartmental models and kinetic systems. In Section 3 the kinetic representation of the studied model class is described, while new results on persistence and periodic behaviour in the time-varying case are proposed in Section 4. Stability analysis results with a family of non-unique logarithmic Lyapunov functions are described in Section 5, and finally, Section 6 summarizes the main results of the paper.

Notations and background
In this section, we describe the basic notations and building blocks of a compartmental system class and chemical reaction networks (CRNs). The notations and overview in this section are based on [27,29].

Compartmental models
Throughout the paper we consider systems containing a set of interconnected compartments and objects (such as ribosomes, particles, molecules, vehicles etc.) moving between them. We assume that the rate of transfer between compartments depends on the amount of objects in the source compartment as well as on the amount of free space in the target compartment. This naturally implies that each compartment has a well-defined finite capacity that limits the amount of modeled quantities that can be contained in the given compartment. We also allow explicit time dependence and in some cases dependence on the amount of objects and free space in other compartments.
For the formal definition, let us consider the set Q = {q 1 , q 2 , . . ., q m } of compartments and the set A � Q × Q of transitions, where (q i , q j ) 2 A represents the transition from compartment q i into q j . Then, the directed graph D = (Q, A) is called the compartmental graph and it describes the structure of the compartmental model. The transitions are assumed to be immediate, thus loop edges are not allowed in the model since they do not introduce additional dynamical terms. Similarly, we do not allow parallel edges between two compartments in the same direction since they can be replaced by a single transition. We say that a (compartmental) graph is strongly connected if there exists a directed path between any two vertices in both directions, and we say that a graph is weakly reversible if it is a collection of isolated strongly connected subgraphs.
For each compartment q i we introduce the sets of donors and receptors, respectively, as that is, the set of donors of a given compartment are the compartments where an incoming transition originates from and the set of receptors are the compartments where an outgoing transition terminates in.

Chemical reaction networks (kinetic systems)
In this subsection we give a brief introduction of kinetic systems based on [18,19], where more details can be found. A chemical reaction network (CRN) contains a set of species S = {X 1 , X 2 , . . ., X N } and the corresponding species vector is given by X = [X 1 X 2 . . . X N ] T . The species of a CRN are transformed into each other through elementary reaction steps of the form where C j ¼ y T j X and C j 0 ¼ y T j 0 X are the source and product complexes, respectively, the vectors y j ; y j 0 2 N N 0 are stoichiometric coefficient vectors and functions K j : R þ N � R þ 7 !R þ are the rate functions with R þ denoting the set of nonnegative real numbers. The matrix Y containing the stoichiometric coefficient vectors as columns is called the stoichiometric matrix. The subspace S � R N spanned by the so-called reaction vectors y j 0 − y j is called the stoichiometric subspace of the CRN. The CRN structure can be uniquely described by a directed graph as follows. For each complex we assign a vertex in the graph and for each elementary reaction step of the form C j ! C j 0 we assign a directed edge between the corresponding vertices. We call the resulting graph the reaction graph of the CRN. The deficiency of the CRN is defined as δ = m − ℓ − s, where m is the number of distinct complexes, ℓ is the number of linkage classes (graph components) in the reaction graph and s is the dimension of the stoichiometric subspace.
Let xðtÞ 2 R þ N denote the state vector of the species as a function of time for t � 0. Based on the above, the dynamics of the CRN is given by We assume that a reaction can only take place if each species of the given reaction have nonzero concentration; that is, we assume that K j ðxðtÞ; tÞ ¼ 0 whenever there exists k 2 supp(y j ) such that x k (t) = 0, where we say that k 2 supp(y j ) if [y j ] k > 0. This property ensures the invariance of the nonnegative orthant (or a part of it). We also presume standard regularity assumptions of the rate functions that guarantee local existence and uniqueness of solutions. Different results in this paper require different sets of such assumptions, thus for the sake of generality they will be specified later. Dynamics of the form of (3) is called persistent if no trajectory that starts in the positive orthant has an omega-limit point on the boundary of R N þ . We note that for any v 2 S ? (where S denotes the stoichiometric subspace) we have that and thus hx, vi is constant. Since v 2 S ? was arbitrary we have that xðtÞ 2 xð0Þ þ S. This shows that the translates of S define invariant linear manifolds for the system. We further define for each p 2 R n þ a positive stoichiometric compatibility class S p ¼ ðp þ SÞ \ R þ n .
A set of ODEs of the form _ x ¼ f ðx; tÞ is called kinetic if it can be written in the form (3) with appropriate rate functions and stoichiometric coefficient vectors.

Kinetic representation
In this section we construct a kinetic representation of the above compartmental system class. To do so, we assign a CRN that incorporates the compartmental structure. This allows the introduction of a system of ODEs of the form (3) describing the time evolution of the compartmental model. Some of the following steps are described in [27] or [29] in a time-invariant setting but here we recall and extend them for convenience.

Kinetic modeling of compartmental transitions
Let us consider a compartmental model D = (Q, A). Let the set of species be S = {N 1 , N 2 , . . ., N m } [ {S 1 , S 2 , . . ., S m } where N i and S i represent the number of particles and available spaces in compartment q i , respectively. To each transition (q i , q j ) 2 A we assign a reaction of the form (see, also [24]) where K ij is the rate function of the transition. Such a reaction represents that during the transition from compartment q i to compartment q j the number of items decreases in q i and increases in q j , while the number of available spaces increases in q i and decreases in q j . Let n i and s i denote the continuous amount of particles and free space in q i , respectively. Based on (3) the dynamics of the system is given by where n and s denote the vectorized form of the variables n i and s i , respectively. It is easy to check that the model class in Eq (6) contains ribosome flow models described in [23] or [15], and extends them in two ways: firstly, the reaction rate function K is not necessarily massaction type and moreover, is time-varying, and secondly, the compartmental graph of the system can be arbitrary (i.e., there can be transitions between any two compartments). Note, that we also allow the transition rates to depend on the amount of objects and free space in other compartments as well, possibly describing inhibitory phenomena. Therefore, we call (6) a generalized time-varying ribosome flow model. Thus, our novel results not only extend the theory of ribosome flow models, but can be applied to other TASEP based transport models [30][31][32][33][34] and other flow models, such as the Traffic Reaction Model of [25] or the Nonlocal Flow Reaction Model of [26]. Finally, we note, that while more complicated network structures may not be biologically relevant in the case of ribosome flows, but can serve as a great tool for the analysis of other flow based physical models, e.g. traffic flows.
Clearly the reaction graph of the assigned CRN of a compartmental model is generally not strongly connected nor weakly reversible even if the compartmental graph is strongly connected. In fact, the reaction graph is weakly reversible if and only if each transition in the compartmental system is reversible. Even though the reaction graph, in some sense, loses the regularities of the compartmental graph, we can explicitly determine its deficiency from the compartmental topology and, as described in [27], CRNs of the form (6) exhibit persistence and stability properties in various senses in the time-invariant case. Proof. For each transition between q i and q j we assign two complexes, namely N i + S j and S i + N j , regardless of the transitions' direction, so reversible reactions do not introduce additional complexes, and thus the number of stoichiometrically distinct complexes is m ¼ 2jÃj. A complex of the form N i + S j is only connected with the complex S i + N j , and thus we have ' ¼ jÃj linkage classes each consisting of exactly two complexes. To find the dimension of the stoichiometric subspace, denoted by s ¼ dim S, observe that the reaction vector of a reaction of the form N i + S j ! N j + S i is

Deficiency of CRNs realizing compartmental models
where e k 2 R 2m denotes the kth unit vector. Again, since y i!j = −y j!i it suffices to consider the undirected graph |D|. Assume that y i!j is such that Then by (7) we have that for each non-zero term of the form c .!l 0 y .!l 0 the right-hand side also contains at least one non-zero term c l 0 !. y l 0 !. , including the terms c i!. y i!. and c .!j y .!j . This shows that the edges corresponding to the reaction vectors of the right-hand side form possibly multiple cycles in |D|. Without the loss of generality we may assume that this subgraph does not contain cycles isolated from (q i , q j ). We have to consider the following cases: 1. First, we assume that the right-hand side is a single chordless cycle and contains the transitions Taking the inner product of unit vectors e i ; e l 1 ; e l 2 ; . . . ; e l r ; e j and yields the system of linear equations: which clearly has one solution where each weight is equal to one. 2. If the right-hand side consists of multiple cycles, then repeatedly using the previous argument we can replace the arcs not containing (q i , q j ) with chords. Note, that if the reaction vector corresponding to the chord is already on the right-hand side, then we just have to modify its coefficient. This method decomposes the right-hand side and will leave us with one chordless cycle containing (q i , q j ), leading back to the previous case with exactly one solution. Repeating the arc substitutions we can see that each arc becomes a chordless cycle with the reintroduced edges and the arising systems of linear equations have exactly one solution.
The first case above shows that the dimension of the stiochiometric subspace reduces by one for each set of reaction vectors that correspond to edges forming a chordless cycle in |D| and the second case shows that is reduced by that exact amount. If σ denotes the number of chordless cycles inQ, then the deficiency of the reaction network can be computed as d ¼ m À ' À s ¼ 2jÃj À jÃj À ðjÃj À sÞ ¼ s.

Linear conservation laws
System (6) exhibits conservation in several senses. First of all, we have that thus the sum of modeled quantities and free spaces in the system is constant along the trajectories of (6); that is, the function H : R 2m 7 ! R defined for x 2 R 2m as is a first integral, where x 1 , x 2 , . . ., x m and x m+1 , x m+2 , . . ., x 2m correspond to the variables n 1 , n 2 , . . ., n m and s 1 , s 2 , . . ., s m , respectively. Our next observation is that _ or after an analogous substitution, as As a consequence of the preceding observations, the functionH : R m 7 ! R, defined for is a first integral for (14), in which case each x i = n i (and similarly for (15) if each x i = s i ). This shows that while the state space of the decomposed systems is ; c m �, for a given initial condition xð0ÞC the trajectories are contained in the (m − 1)-dimensional manifold (hyperplane) defined by fx 2CjHðxÞ ÀHðxð0ÞÞ ¼ 0g: For a generalized ribosome flow define c ¼ P n i¼1 c i and for r 2 [0, c] let L r �C be the level set of H corresponding to r; that is, Example 1.1: (generalized) RFMR. As a small example let us consider a Ribosome Flow Model on a Ring (RFMR) [13] with three sites. The underlying compartmental model is given The topology is shown in Fig 1. The corresponding CRN has the following species and reactions: It is easy to see that, indeed, the reaction graph is not weakly reversible and its deficiency is one. The dynamics of the model in the full state space is given by (6)  which can be rewritten in the reduced state space based on (14) as _ n 1 ¼ K 31 ðn; c ðmÞ À n; tÞ À K 12 ðn; c ðmÞ À n; tÞ _ n 2 ¼ K 12 ðn; c ðmÞ À n; tÞ À K 23 ðn; c ðmÞ À n; tÞ _ n 3 ¼ K 23 ðn; c ðmÞ À n; tÞ À K 31 ðn; c ðmÞ À n; tÞ: In a classical RFMR each c i = 1 and each transition-rate K ij depends only on n i and c i − n i = 1 − n i , and follows the mass-action law. In an RFMR with different site sizes (RFMRD) [15] we allow arbitrary site sizes, in which case the above equation can be written as

Analysis of persistence and stability
In this section we show that systems of the form (6) exhibit various interesting dynamical properties that can be characterized under different assumptions of the transition rate functions. First we will consider time-invariant systems to demonstrate the regularity of equilibria. Then we return to time-varying systems to generalize the results of [27].

Equilibria of time-invariant systems
In this subsection we assume that the K ij ðn; s; tÞ rate functions are continuously differentiable and only depend on the variables n i and s j in a nondecreasing manner; that is, we assume that K ij ðn; s; tÞ � K ij ðn i ; s j Þ for each i and j. Then the results [27] show that a system of the form (14) is cooperative (the name also highlights the importance of the exclusion of inhibitory phenomena), is (strongly) monotone and each level set L r contains a unique globally (relative to its level set) asymptotically stable steady state. This implies that the steady states form a linearly ordered set. For i = 1, 2, . . ., m let e i : [0, c] 7 ! [0, c i ] denote the ith coordinate function of the steady state; that is, let where n(0) 2 L r is arbitrary and ρ(t, n(0)) denotes the solution at time t with ρ(0, n(0)) = n(0).
Clearly each e i is continuous and the monotonicity of the system also shows that each e i function is strictly increasing; that is, they are differentiable almost everywhere and their derivative are positive. Example 1.2: Equilibria of generalized RFMR. Let us consider a generalized version of the RFMR in Fig 1. Its time evolution in the reduced state space is given in the form (14) as and for the simulations we set capacities c 1 = 5, c 2 = 25, c 3 = 50. The rate functions in the different cases are assumed to have the form K ij ðn i ; c j À n j Þ ¼ k ij n i ðc j À n j Þ (corresponding to the classical RFMRD) or to be rational functions of the form for some l > 0 with k 12 = 100, k 23 = 40, k 31 = 60. Fig 2 shows the equilibrium curves for these rate functions with various l values. Example 2: Not strongly connected model. Let us consider consider a not strongly connected compartmental model given by The topology is shown in  The corresponding CRN has the following species and reactions: The dynamics of the system in the reduced state space is given by Since the compartmental graph is not strongly connected the persistence and stability results of [27] are not applicable. However, empirical results show that the long-time behaviour of the system still exhibits some regularity, which can be divided into two cases base on the initial values of the system: and n 1 (t) and n 2 (t) will converge to the unique equilibrium on the level set Note that since D 0 is strongly connected, the results of [27] and the above investigation can be applied.
For the simulations we set c 1 = c 2 = c 3 = 100. The rate functions in the different cases are assumed to have form K ij ðn i ; c j À n j Þ ¼ k ij n i ðc j À n j Þ (corresponding to mass-action kinetics) or to be rational functions of the form Fig 4 shows the equilibrium curves for these rate functions with various l values. As described by the above cases we see that until the sum of the initial value exceed the capacity of the q 1 compartment the equilibrium lies on the n 1 axis. After that the equilibrium lies on the plane fn 1 ¼ c 1 g � R 3 and since D 0 is strongly connected we have that the coordinate functions of the equilibria e 2 (r) and e 3 (r), restricted to the set [c 1 , c], are continuous and strictly increasing. We note that while the system is not strongly connected it exhibits many similar qualitative properties as strongly connected models. For example, for initial values satisfying H(n(0)) > c 1 the system is Lyapunov stable as described in [29].   Based on the initial value and the exact compartmental structure the following phenomena can be observed: • Traps (and only traps) can become full, thus possibly creating new traps.
• Sources (and only sources) can become empty, thus possibly creating new sources.
• After a sufficient number of traps are filled and sources are emptied, the compartmental graph D is decomposed into isolated strongly connected components; that is, the resulting graph is weakly reversible, in which case the results of [27] can be applied.
While these observations are elementary and show that the system is stable, the equilibria are clearly non-unique with respect to the total mass of the network and in general it is not straightforward to predict from the initial value which components will fill and empty.

Persistence
In this subsection we consider time-varying generalized ribosome flows of the form (6) only under mild regularity assumptions described by the following theorem, which is based on the results of [35] but the statements are rephrased to be more aligned with our framework. For the definition of notions related to Petri nets (e.g. siphons) and their exact connection with CRNs we refer to [27,35]. To verify condition (i) we would, in general, need to enumerate all siphons of the CRN, which is well-known to be an NP-hard problem. However, in our recent paper [27] we explicitly characterized the siphons of a CRN assigned to a strongly connected compartmental models in the time-invariant case. However, one can observe that conditions (i) and (ii) of 4.2 are independent of the choice of transition rates and even independent from whether the system is time-invariant or not; that is, our results, formulated in the following theorem, hold for timevarying compartmental systems as well. Then the conclusions of Section 3.3 show that conditions (i) and (ii) are satisfied by virtue of the first integrals (16) and (13), respectively.
It is not straightforward to determine exactly what types of reaction rates satisfy condition (iii). For the sake of specificity, we characterize a class of reaction rates of special interest which can be written in the following form where we assume that the transformations y i ; n j 2 C 1 ðRÞ are nondecreasing, have θ i (0) = ν j (0) = 0 and satisfy R 1 0 jlog y i ðrÞjdr < 1 and R 1 0 jlog n j ðrÞjdr < 1 for each i, j = 1, 2, . . ., m. We also assume that the functions C ij take the form where r ð1Þ ; r ð2Þ 2 N m and a r ð1Þ ;r ð2Þ 2 R þ . We further assume that for k ij (t) there exist k ij ; k ij > 0 such that k ij � k ij ðtÞ � k ij for all t � 0. In this case we have which are clearly monotonous in the sense of Theorem 4.2, and thus condition (i) is satisfied and the system is persistent. The denominator of (29) contains positive terms which can be interpreted as the inhibitory effect of other species, and the time-varying coefficient k ij (t) introduces the dependence of the system parameters on various factors such as temperature or the dynamical behaviour of other species that are not explicitly modeled as state variables. This class of rate functions contains many well-known examples, demonstrating the range and flexibility of reaction rates of the above form: 1. Setting each θ i (n i ) = n i and ν j (s j ) = s j and C ij (n, s) = 0 we obtain the case of classical massaction kinetics with time-varying rate coefficients: K ij ðn; s; tÞ ¼ k ij ðtÞn i s j .
2. Setting each θ i (n i ) = n i and ν j (s j ) = s j and C ij (n, s) = l 2 − 1 + ln i + ls j + n i s j for some l > 0 yields K ij ðn; s; tÞ ¼ k ij ðtÞ n i s j ðl þ n i Þðl þ s j Þ ð32Þ corresponding to simple saturating kinetics described by the Monod equation.
3. The previous example can also be obtained by setting y i ðn i Þ ¼ n i lþn i and n j ðs j Þ ¼ s j lþs j and C ij (n, s) = 0, showing that (29) is not unique. Notice however, that for fixed θ i , ν j transformations the function C ij , and thus the fraction itself, is unique. For this example we set c 1 = c 2 = c 3 = 100, l = 100 and K 12 ðn 1 ; c 2 À n 2 ; tÞ ¼ k 12 ðtÞ n 1 ðc 2 À n 2 Þ ðl þ n 1 Þðl þ c 2 À n 2 Þ ; K 23 ðn 2 ; c 3 À n 3 ; tÞ ¼ k 23 ðtÞ n 2 ðc 3 À n 3 Þ ðl þ n 2 Þðl þ c 3 À n 3 Þ ; where the coefficient functions are considered to be exponentially decaying perturbations of the nominal values of the form As a comparison let us consider the solutionñðtÞ of the time-invariant system with the above nominal values. We can observe that since the time dependent terms are exponentially decaying and both systems evolve on the same linear manifold, the systems tend to the same equilibrium, as expected.

Stability of the solutions for periodic transition rates
In this section we investigate the periodic behaviour of the generalized ribosome flows based on the ideas of [36]. Let us consider a generalized ribosome flow in the reduced state space of the form (14) with transition rates of the form (29) and assume that the transition functions are C 1 and periodic with the same period (but having possibly different phases). Write (14) as _ n ¼ Fðt; nÞ and assume that the right-hand side satisfies the following monotonicity condition: F i (t, x) � F i (t, y) for any two distinct points x; y 2C such that x i = y i and x j � y j for j 6 ¼ i. This condition is satisfied if, for example, the transition rates are such that C ij � 0; that is, if there are no inhibitory phenomena. Then the system phase locks (or entrains) with the periodic excitations.

Lyapunov stability analysis
In this section we show that generalized ribosome flows with reaction rate functions of the form (29) with piecewise locally Lipschitz k ij (t) coefficients satisfy a certain notion of robustness to the changes in the time-varying rate functions that can be traced back to the input-tostate stability of rate-controlled biochemical networks thoroughly investigated in [22]. The main difficulty in applying these results lies in the aforementioned fact that the CRN assigned to a compartmental model is generally not weakly reversible and its deficiency is generally not zero (see, Theorem 3.1) even if the compartmental topology is strongly connected. In order to circumvent this, we will perform a model reduction and rewrite (14) by factoring out appropriate terms. Let us first recall the most important notions and results of [22]. Consider the system corresponding to a CRN with R reactions where the nonnegative functions u ij are piecewise locally Lipschitz with a finite number of discontinuities and the stoichiometric coefficient vectors y i , y j are as described in 2.2. Motivated by control designs for ribosome flow models [39] we introduce such time dependence not only to handle some uncertainty originating from fluctuating external factors but to measure the robustness of the system to certain control inputs. In this section, however, we restrict the conditions on the transformation functions y i : R þ 7 ! ½0; 1Þ. Namely, we assume that (a) θ i is real analytic, R 1 0 jlog y i ðrÞjdr < 1 (d) θ i is strictly increasing and onto the set [0, σ i ) for some σ i 2 [0, 1), (e) lim t!log s i R t a r À 1 i ðrÞdr À pt ¼ 1 for any a < log σ i and any constant p > 0, where ρ i = log θ i .
Before continuing with the definitions, we consider the case when u(t) is a constant matrix A. We assume that A has nonnegative entries and is irreducible; that is, the underlying reaction graph is strongly connected. We denote the set of such A matrices as A. Then the equilibria of _ x ¼ f ðx; AÞ can be divided into the sets of boundary equilibria and positive equilibria: Then, the result [22, Theorem 2.1] (and also [40,Theorem 2]) shows that if there are no boundary equilibria in any positive class, then each positive class contains a unique globally (relative to the positive class) asymptotically stable positive equilibrium. Denote the unique positive equilibrium in the same class as x 0 as xðx 0 ; AÞ and notice that Definition 5. 1 We define the following function classes: (i) A function a : R þ 7 ! R þ is said to be of class K if it is continuous, strictly increasing and has α(0) = 0.
(ii) The subset of unbounded functions of class K are denoted by K 1 . 0 and β(r,.) is strictly decreasing to zero for all r > 0.
We consider nonnegative time-varying inputs such that at any time instant the reaction graph is strongly connected; that is, the input-value set U is a subset of A. Furthermore, let k.k 2 denote the spectral norm induced by the Euclidian norm and for u : Definition 5.2. A system _ x ¼ f ðx; uÞ is uniformly input-to-state stable (ISS) with input-value set U if for every compact set P � E and every compact set F � R þ n containing P, there exist functions β = β P of class KL and ϕ = ϕ P of class K 1 such that, for every x o 2 P \ E u 0 ;þ for some u 0 2 U we have that holds for each u : R þ 7 ! U input and every initial condition According to the above definition we say that a system is ISS if it is globally asymptotically stable in the absence of external inputs and if its trajectories are bounded by an appropriate function of the input. In some sense this definition is intended to capture the idea of "bounded input bounded output" stability, since for bounded u input (u − u 0 to be more precise) the trajectories will remain in a ball and, in fact, approach the ball �ðku À u 0 k U Þ as t increases [41].
We assume that there exists a uniform lower bound on the parameters; that is, we consider input-value sets of the form We also recall that the input functions are piecewise locally Lipschitz in time with a finite number of discontinuities, thus we introduce Then the main Theorem of [22] states: Theorem 5.3. Consider the system (38) with and suppose that is is mass-conservative; that is, there exists v 2 R n þ such that v T f(x, u) = 0 for all x 2 R þ n and u 2 A. Then the system with input maps u 2 W is uniformly ISS with input-value set U 2 .
The proof relies on the candidate ISS-Lyapunov function (for the definition of which and for the exact connection with ISS stability we refer to [22]) which, for mass-action systems, yields the classical entropy-like Lyapunov function wellknown from the theory of chemical reaction networks, see (65). We note that Vðx; xÞ is uniquely determined by the θ i functions and does not depend explicitly on the reaction/compartmental structure or the time-varying u ij (t) functions; that is, it is universal in the sense of [42]. [29] for more details), in which case the above Lyapunov function can be applied.

Factorization of the transition rates
Let us consider a generalized ribosome flow in the reduced state space of the form (14), in this case given by K ji ðn; c À n; tÞ À X j2R i K ij ðn; c À n; tÞ ¼ X j2D i k ji ðtÞ y j ðn j Þn i ðc i À n i Þ 1 þ C ji ðn; c ðmÞ À nÞ À X j2R i k ij ðtÞ y i ðn i Þn j ðc j À n j Þ 1 þ C ij ðn; c ðmÞ À nÞ : Notice that we can naturally factor some terms of the transition rates into the time-varying coefficient as Then (46) can be rewritten as This equation can be clearly embedded into the class of strongly connected systems of the form (38), since the reaction graph of (48) consists of species S = {N 1 , N 2 , . . ., N m }, has the m × m identity matrix as its stoichiometric matrix and for each transition (q i , q j ) 2 A we assign a reaction of the form and thus the system of differential equations can be written as where the elements ofÃ k are given by Note that the fractions n j ðc j À n j Þ 1þC ij ðn;c ðmÞ À nÞ are differentiable (and thus Lipschitz) and each k ij (t) is piecewise locally Lipschitz, hence eachk ij ðtÞ is piecewise locally Lipschitz. This shows that generalized ribosome flows can be embedded into the class of rate-controlled biochemical networks described in [22] in a way that preserves the compartmental structure; that is, the reaction graph of (50) is topologically isomorph to the compartmental graph. In particular if the compartmental model is strongly connected, then the reaction graph of the reduced system (50) is strongly connected as well. Furthermore, combining the persistence of the system with Remark 4.4 we find thatÃ k 2 W, and thus Theorem 5.3 ensures input-to-state stability.

Quasi-LTV factorization
A classical argument shows that the model reduction above can result in a Linear Time-Varying (LTV) system [6]. Consider an FðxÞ 2 C k ðRÞ nonnegative function such that F(0) = 0, where k � 1. Then for the function F(rx) we have and thus and since F(0) = 0, we find that F(x) = xf(x). Note, that the calculation also shows that f 2 C kÀ 1 ðRÞ. Since θ i is real analytic we have that y i ðn i Þ ¼ŷ i ðn i Þn i for someŷ i real analytic function. Then (48) can be rewritten as wherek ij ðtÞ ¼ Similarly as before, the reaction graph of (54) consists of species S = {N 1 , N 2 , . . ., N m }, has the m × m identity matrix as its stoichiometric matrix and for each transition (q i , q j ) 2 A we assign a reaction of the form and thus the system of differential equations can be written as where the elements ofÂ k are given by Again, eachk ij ðtÞ is piecewise locally Lipschitz, thus for strongly connected compartmental models Theorem 5.3 ensures input-to-state stability via Remark 4.4.
Then the CRN corresponding to (64) has the following species and reactions: which is strongly connected and isomorph to the compartmental model in Fig 1. We arrive at the same conclusion if we instead use the functionŝ to rewrite (59) as Note that the quasi-LTV factorization might be more complicated in some cases, but the construction described in Section 5.2 guarantees its existence.

Induced family of Lyapunov functions
The above investigation demonstrates that generalized ribosome flows can be embedded into rate-controlled biochemical networks in at least two different ways, where each embedding induces a different Lyapunov function of the form (45). Thus, in general, we may use at least two different Lyapunov functions governing the same dynamics. To characterize their exact relation, consider a factored system of the form (50) with its ISS-Lyapunov function Vðn; nÞ. The quasi-LTV representation of the system admits an ISS-Lyapunov function of the form  [43]. While in general we are restricted to the above factorizations, in some special cases we may use a whole family of factorizations and corresponding Lyapunov functions. To illustrate this, consider an example when each y i ðrÞ ¼ r a i ðlþrÞ b i for some l > 0 and a i 2 N, b i 2 N 0 , a i � b i (these properties ensure that the functions θ i are nondecreasing). Then, after the factorization described in Section 5.1, the Lyapunov function (45) becomes V ðl;a;bÞ ðn; nÞ ¼ We emphasize that (45) only depends on the θ i functions, in this case parametrized with the l, a i , b i values; that is, it is independent of the network structure and transition rate coefficients.
We can also perform the factorization y i ðrÞ ¼ỹ i ðrÞ râi ðlþrÞb i withâ i 2 N,â i < a i ,b i 2 N 0 ,â i �b i yielding the Lyapunov function V ðl;â;bÞ of the same form as in (67). This shows that the parameters a and b can be freely (apart from the constraints above) chosen in (67). We may also observe some interesting behaviour at the extrema of the parametersb and l, namely, that if we choose eachb i ¼ 0 then the Lyapunov function in (67) is independent of l; that is, we have that Moreover, letting l!1 yields the convergence where V LTV i is defined in (65). Example 1.5: Family of Lyapunov functions of a generalized RFMR. Let us again consider a generalized version of the RFMR in the reduced state space from Fig 1. For a given initial condition n 0 we can substitute n 3 = H(n 0 ) − n 1 − n 2 , and thus the Lyapunov function restricted to the manifold fHðnÞ ¼ Hðn 0 Þg can be seen as a two dimensional function with local coordinates n 1 and n 2 .
We set the capacities as c 1 = c 2 = c 3 = 100 and k 12 = 100, k 23 = 60, k 31 = 20. The system has transition rates as described above with each a i = b i = 3; that is, we have that  � � Example 4: Competition for ribosomes in the cell. In this example we introduce a set of generalized ribosome flows connected by a finite pool of ribosomes to model competition in the cell. We follow [14], where the authors introduced a model for simultaneous translation and [16], where the authors generalized the model to include premature drop-off and attachment effects modeled with Langmuir kinetics. We will focus on the latter case and show that with a slight modification it can be formalized as a generalized ribosome flow model with a clear and natural compartmental interpretation. This demonstrates the usefulness and modeling power of generalized ribosome flows as one can prove various properties of many existing models of different conceptual levels. Moreover, our results show that many qualitative properties of the system carry over to more general settings, e.g. when the translation, drop-off and attachment rates are modeled with more sophisticated functions or when some (or all) rates are time-dependent.
For the sake of simplicity we will present this example in the reduced state space. Let us consider N mRNAs consisting of m 1 , m 2 , . . ., m N number of sites. Let n j i denote the continuous amount of ribosomes in the ith site of the j mRNA stand and let c j i denote its capacity. Let c z denote the capacity of the pool and n z denote the amount of ribosomes in the pool. For the sake of notational simplicity let n j 0 and n j m j þ1 also denote n z and similarly for the capacities. Let the translation rate functions from the ith site the to (i + 1)th site on the jth mRNA be denoted as K j iðiþ1Þ . Finally, let the detachment and attachment rates at the ith site of the jth mRNA be denoted respectively as K j iz and K j zi . The attachment rate to the first site and the detachment rate from the last site will be called initiation rate and production rate, respectively. Then the dynamics of the model is given by: ; tÞ þK j zi ðn z ; c j i À n j i ; tÞ À K j iz ðn j i ; c z À n z ; tÞ; Thus, indeed, simultaneous translation with a finite pool can be described by a generalized ribosome flow. Clearly the following function defines a linear first integral and is a crucial factor in the dynamical analysis of the system. Remark 5.7. In [16] the authors consider the following special case: • the capacity of each site is one; that is, each c j i ¼ 1, • the translation rate are time-invariant and obey the mass-action law; that is, each K j iðiþ1Þ ðn j i ; 1 À n j iþ1 ; tÞ ¼ l j i n j i ð1 À n j iþ1 Þ for some l j i > 0, • the initiation and attachment rates are time-invariant and are given by K j zi ðn z ; 1 À n j i ; tÞ ¼ b j i G j ðzÞð1 À n j i Þ for some b j i � 0 and G j (z) continuously differentiable strictly increasing function with G j (0) = 0, • the drop-off rates are time-invariant and are given by K j iz ðn j i ; c z À n z ; tÞ ¼ a j i n j i for some a j i � 0.
Since the drop-off rates are donor controlled the pool does not have a predefined capacity and the amount of ribosomes in the pool are only bounded by H(n(0)). Therefore, this special case does not fit in our compartmental framework (although, as most of our results are a consequence of the linear first integral combined with the cooperativity of the system they can be generalized to include donor controlled terms as well). It is assumed that the authors consider this case to capture the fact that the capacity of the pool might be several orders higher than the actual number of ribosomes, and thus the dependence on the available space in the pool may be negligible. However, some physical meaning is lost with this assumption and it might in fact lead to less precise simulations.
To see this, let us consider a network with N = 10 mRNAs with m = 5 sites. For the sake of simplicity let l j i ¼ b j 1 ¼ a j 5 ¼ 1 for each i and j, and assume that there are no premature drop-offs and attachments. We consider initation rates G j (z) = z, G j (z) = tanh(z) and G j (z) = z 2 and set c z = 10 4 . Since the equilibrium is unique on the level sets of the first integral we set each n j i ¼ 0 and we only change n z (0). Fig 9 shows the ratio of the steady state of the pool in the case of donor controlled and mass-action production rates as we increase the ratio n z ð0Þ c z from 5 � 10 −2 to 1. As expected, the steady state ratio is close to one for saturating rate functions and for n z (0) � c z . However, the ratio can get higher when the total number of ribosomes have the same magnitude as the capacity; that is, the inaccuracy of the donor controlled kinetics increases. While this assumption might be valid for realistic parameters of ribosome flows in other TASEP based flow models (especially with non-saturating kinetics) it might be crucial to model these transitions accurately.
Effect of the total number of ribosomes. In the next simulation we follow [16] and we consider a single mRNA strand with m 1 = 3 sites. The initiation rate is set to b 1 1 ¼ 1 while the attachment rates are b 1 2 ¼ 0:1 and b 1 3 ¼ 0. The drop-off rates and production rate are set as a 1 1 ¼ 0, a 1 2 ¼ 0:01, a 1 3 ¼ 1. We assume that the translation rates obey the mass-action law with each l 1 i ¼ 1. We set the initial values to n 1 j ¼ 0 and n 0 (z) = c z as before. Fig 10 shows the steady state of the system as we increase c z from 0 to 5 for various rate functions. One can see that in each case the mRNA saturates as we increase the number of ribosomes and the rest of the ribosomes are accumulated in the pool. Finally, the same effect as in Fig 9 can be observed; that is, the donor controlled detachment rates shift the steady state of the pool to higher values.
We again emphasize the versatility of generalized ribosome flows as the initiation, translation, production, attachment and detachment rate function can be different on each site. For example let us consider a particular mRNA strand with saturating initation and attachment  rates given by K 1 zi ðn z ; n 1 i Þ ¼ b 1 i tanhðn z Þðc 1 i À n 1 i Þ, with mass-action translation rates and with production and drop-off rates given by K 1 iz ðn 1 i ; n z Þ ¼ a 1 i � n 1 i 1þn 1 i � n 3 z . Fig 11 shows evolution of the steady states as we increase n z (0) = c z as before. As expected the steady states of the mRNA sites are moved to lower values.

Conclusions
The dynamical properties of generalized ribosome flow models models with arbitrary compartmental graph structure and general time-varying transition rates were studied in this paper. The analysis is based on the deterministic CRN representation of such systems which has a transparent physical meaning by tracking the amounts (concentrations) of available objects and free spaces, respectively, in each compartment. Our framework includes several important models from the literature including the RFMR [13], and its generalizations like the RFMRD [15]. As demonstrated in Example 4, our framework can describe complex phenomena like competition for ribosomes in a cell through a set of tubular flow models connected with a pool of finite capacity. The obtained model (with a slight modification due to physical considerations) includes previously published pool models such as [16]. It was shown that the deficiency of the obtained kinetic model form is equal to the number of chordless cycles in the undirected reaction graph of the system. Furthermore, it was proved that time-varying generalized ribosome flows are persistent under mild regularity assumptions on the transition rates, and a wide set of reaction rates satisfying this assumption was characterized, containing well- known examples such as mass-action type rates. It was shown that the studied models can be embedded in at least two ways into the class of rate-controlled biochemical networks originally described in [22]. This embedding allows us to prove stability with entropy-like logarithmic Lyapunov functions known from the theory of CRNs. It was illustrated that the non-unique factorization of the rate functions gives rise to a whole family of various possible Lyapunov functions. Finally, periodic model behaviour was also studied, where we showed that trajectories with the same overall initial mass and periodic transition rates having the same period (but possibly different phase) converge to a unique periodic solution. The generality of our setting allows the efficient extension of the proposed results to ribosome flows open to the environment, for example as in [44,45]. It also allows the generalization of other qualitative results in the literature, for example translation rate maximization [15,46]. We emphasize that the persistence and stability results, and even the basic logarithmic Lyapunov function candidate in Eq (45) are independent of the potentially uncertain time varying parameters (rate coefficients) of the model. Besides the theoretical aspects, these improvements may efficiently support structural design or control of compartmental models with bounded capacities, which is planned to be addressed in our future work.